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1.  Summary  of  research 

We  have  performed  research  on  the  algorithm  development,  analysis,  implemen¬ 
tation  and  application  for  high  order  weighted  essentially  non-oscillatory  (WENO)  fi¬ 
nite  difference  and  finite  volume  schemes,  discontinuous  Galerkin  (DG)  finite  element 
method,  and  related  methods,  for  solving  computational  fluid  dynamics  (CFD)  prob¬ 
lems  and  other  applications  containing  strong  shock  waves  and  complicated  smooth 
region  structures.  The  goal  of  this  research  is  to  help  obtaining  more  robust,  cost 
effective,  and  reliable  numerical  tools  for  computational  fluid  dynamics  problems  and 
other  applications,  that  are  of  interest  to  AFOSR. 

There  are  17  publications  in  refereed  journals,  which  have  appeared  or  have  been 
accepted  for  publications  and  have  quoted  partial  support  by  this  grant,  see  section 
3. 

In  [1]  (the  numbers  here  and  below  refer  to  that  of  publications  listed  in  section 
3),  we  have  explored  an  alternative  formulation  of  finite  difference  WENO  schemes 
with  Lax-Wendroff  time  discretization  for  solving  hyperbolic  conservation  laws.  Even 
though  this  alternative  formulation  is  more  expensive  than  the  standard  formulation, 
it  does  have  several  advantages.  The  first  advantage  is  that  arbitrary  monotone  fluxes 
can  be  used  in  this  framework,  while  the  traditional  practice  of  reconstructing  flux 
functions  can  be  applied  only  to  smooth  flux  splitting.  This  allows  us  to  choose 
monotone  fluxes  or  approximate  Ricmann  solvers  with  smaller  numerical  viscosity. 
The  second  advantage  is  that  a  narrower  effective  stencil  is  used  compared  with  pre¬ 
vious  high  order  finite  difference  WENO  schemes  based  on  the  reconstruction  of  flux 
functions,  with  a  Lax-Wendroff  time  discretization.  The  third  and  major  advantage  is 
that  free-stream-preserving  high  order  finite  difference  WENO  scheme  can  be  designed 
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via  this  formulation,  which  keeps  conservation,  free-stream-preserving,  and  high  order 
accuracy  simultaneously.  This  advantage  is  fully  explored  with  numerical  experiments 
in  [3].  The  scheme  in  [3]  has  attracted  a  lot  of  attention  from  CFD  applications,  as  it 
is  a  well-known  difficult  task  to  have  conservation,  free-stream-preserving,  and  high 
order  accuracy  fulfilled  at  the  same  time  for  high  order  nonlinear  finite  difference 
schemes. 

In  [2],  a  flux  correction  approach  is  used  to  design  high  order  conservative  finite 
difference  and  finite  volume  schemes  which  are  positivity-preserving  for  density  and 
pressure  in  Euler  systems.  Such  methods  are  proved  to  maintain  high  order  accuracy 
under  mild  assumptions.  Computational  results  are  provided  in  [2]  to  demonstrate 
the  good  performance  of  these  schemes. 

In  [4],  we  study  spectral  collocation  methods  for  functions  which  are  analytic  in  the 
open  interval  but  have  singularities  at  end-points.  Naive  summation  of  the  spectral 
partial  sum  leads  to  divergence  near  the  end  points  and  also  poor  convergence  rate 
inside  the  interval.  Previous  techniques  in  post-precessing  to  overcome  this  difficulty 
do  not  apply  because  of  the  end-point  singularities.  We  propose  new  analysis  with 
suitable  post-precessing,  which  leads  to  exponential  accuracy  in  the  maximum  norm 
everywhere  in  the  interval  using  only  the  information  of  point  values  at  standard 
collocation  points  of  such  functions.  The  proof  is  constructive  and  provides  guidance 
for  practical  implementation  of  the  post-processor.  Numerical  experiments  are  given 
verifying  the  spectral  convergence  after  post-processing.  In  [8],  this  technique  and 
its  analysis  are  generalized  to  Galerkin  approximation  and  to  the  solution  of  linear 
PDEs.  The  techniques  in  [4]  and  [8]  have  already  been  used  by  other  scientists  to 
deal  with  interfaces  in  kinetic  problems. 

In  [5],  we  design  first  order  and  higher  order  Lagrangian  schemes  for  multi- material 
compressible  flows,  which  are  positivity-preserving  for  density  and  internal  energy. 
Lagrangian  schemes  are  particularly  useful  for  solving  multi-material  compressible 
flows,  as  they  keep  the  sharp  interface  between  different  materials  during  the  evolu¬ 
tion.  However,  stability  for  Lagrangian  schemes  is  much  more  difficult  to  maintain 
than  for  Eulerian  schemes.  Especially  for  flows  with  general  equations  of  state,  even 
first  order  positivity-preserving  schemes  in  multi-dimensions  were  not  available  before. 
We  have  constructed  not  only  first  order  positivity- preserving  Lagrangian  schemes  in 
multi-dimensions,  but  also  higher  order  finite  volume  schemes  with  the  same  property. 
Numerical  results  demonstrate  the  good  performance  of  these  schemes  for  extreme 
flows,  which  lead  to  blowups  of  the  scheme  without  the  positivity-preserving  tech¬ 
niques. 

Recently,  a  class  of  schemes  termed  correction  procedure  via  reconstruction  (CPR) 
has  been  designed  and  used  by  aerospace  engineers,  who  found  the  schemes  perform¬ 
ing  well  on  both  structured  and  unstructured  meshes.  However,  for  simulations  with 
strong  shocks,  numerical  oscillations  may  lead  to  nonlinear  instability.  In  [6],  we 
have  adapted  a  simple  WENO  limiter  which  uses  information  only  from  immediate 
neighboring  cells,  originally  designed  for  DG  schemes,  to  the  CPR  schemes  on  un- 
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structured  meshes.  Numerical  experiments  indicate  the  good  performance  of  this 
simple  limiter  which  both  keeps  high  order  accuracy  in  smooth  regions  and  controls 
spurious  oscillations  near  shocks.  Because  the  limiter  uses  information  only  from 
immediate  neighboring  cells,  it  does  not  destroy  the  original  local  data  structure  of 
the  DG  scheme  and  hence  it  maintains  the  major  advantage  of  DG  schemes  in  high 
parallel  efficiency. 

In  [9],  a  well-designed  adaptive  mesh  refinement  (AMR)  method  based  on  high 
order  finite  difference  WENO  schemes  is  developed  and  applied  to  simulate  three- 
dimensional  detonation  waves.  Such  simulations  require  large  computational  re¬ 
sources  because  of  the  high  resolution  requirement  for  the  detailed  detonation  struc¬ 
ture.  The  efficient  parallel  AMR- WENO  method  provides  a  good  tool  for  these 
detonation  simulations. 

In  [10],  a  predictive  dynamic  user-optimal  traffic  flow  model  in  the  continuum 
space  with  anisotropic  condition  is  developed,  which  involves  a  coupled  system  of 
a  forward-in-time  conservation  law  for  the  traffic  density  and  a  backward-in-time 
Hamilton- Jacobi  equation  for  the  cost  potential.  This  coupled  system  is  very  inter¬ 
esting  mathematically,  in  terms  of  suitable  assumptions  for  its  well-posedness.  More 
importantly,  it  provides  a  big  challenge  for  its  numerical  simulation,  as  an  iterative 
process  would  be  needed  (alternatively  from  forward  in  time  with  the  Hamilton-Jacobi 
solution  frozen,  followed  by  backward  in  time  with  the  conservation  law  solution 
frozen).  We  have  explored  such  an  iterative  numerical  method  in  [10]  and  tested  it 
for  this  complex  system. 

The  so-called  inverse  Lax-Wendroff  procedure,  which  is  a  technique  to  obtain  sta¬ 
ble  and  accurate  numerical  boundary  conditions  when  a  cartesian  mesh  is  used  to 
solve  a  complex  geometry  problem  (hence  the  grids  are  not  aligned  with  the  compu¬ 
tational  domain  boundary),  was  proposed  by  the  PI  with  his  students  a  few  years 
ago.  This  method  has  now  found  a  wide  usage  in  applications.  In  [11],  we  give  a 
stability  analysis,  using  both  the  GKS  (Gustafsson,  Kreiss  and  Sundstrom)  theory 
and  an  eigenspectrum  analysis,  for  this  procedure  applied  to  upwind-biased  finite 
difference  schemes  for  hyperbolic  conservation  laws.  This  study  provides  a  theoret¬ 
ical  guidance  to  determine  the  minimum  number  of  terms  needed  in  this  procedure 
to  guarantee  stability,  thus  allowing  us  to  design  the  most  efficient  yet  guaranteed 
stable  boundary  treatments  for  such  equations.  In  [17],  similar  stability  analysis  is 
performed  for  the  central  schemes  approximating  the  diffusion  equations.  In  [14], 
we  have  extended  this  technique  to  handle  convection- diffusion  equations  including 
compressible  Navier-Stokes  equations,  and  have  designed  a  procedure  which  works 
well  across  the  convection-dominated  to  the  diffusion-dominated  regimes.  This  is  a 
highly  non-trivial  task  as  boundary  conditions  play  very  different  role  for  diffusion- 
dominated  regimes  and  convection-dominated  regimes,  while  our  numerical  proce¬ 
dure  is  expected  to  automatically  adjust  itself  to  yield  stable  and  accurate  algorithms 
across  these  two  regimes.  The  inverse  Lax-Wendroff  procedure  has  already  been  used 
in  several  applications  and  is  expected  to  play  more  important  roles  in  the  future. 
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Further  development  of  this  technique  is  now  underway. 

In  [12],  we  have  obtained  a  bound- preserving  DG  method  for  solving  the  relativis¬ 
tic  hydrodynamics,  which  guarantees  positive  density  and  pressure  and  the  upper 
bound  of  flow  speed  by  the  speed  of  light.  This  conservative  scheme  is  then  L 1  stable 
for  all  the  conserved  variables,  which  is  rare  for  high  order  schemes  solving  hyperbolic 
systems. 

A  novel  technique  to  use  the  idea  of  DG  schemes  to  design  a  stable,  conservative 
and  high  order  accurate  algorithm  to  solve  hyperbolic  conservation  laws  on  an  ar¬ 
bitrarily  distributed  point  clouds  is  designed,  analyzed  and  tested  in  [15].  By  using 
the  Voronoi  technique  and  by  introducing  a  grouping  algorithm,  we  divide  the  com¬ 
putational  domain  into  non-overlapping  cells.  Each  cell  is  a  polygon  and  contains  a 
minimum  number  of  the  given  points  to  ensure  accuracy.  We  carefully  select  points 
in  each  cell  during  the  grouping  procedure,  and  hence  are  able  to  interpolate  or  fit 
the  discrete  initial  values  with  piecewise  polynomials.  The  DG  framework  allows  us 
to  design  stable  and  accurate  evolution  algorithm  on  such  polygonal  cells,  which  have 
different  number  of  sides  for  different  cells  and  the  cells  may  not  be  convex,  causing 
significant  difficulty  for  traditional  finite  element  methods.  Our  algorithm  provides 
good  numerical  results  on  the  preliminary  test  cases  in  [15].  Further  development  of 
this  technique  is  now  underway. 

In  [16],  a  high  order  positivity-preserving  DG  method  is  designed  for  the  radia¬ 
tive  transfer  equations,  which  involves  a  completely  different  strategy  from  previous 
techniques  as  the  algorithm  is  a  spatially  marching  type  for  steady  states  or  implicit 
type  for  time-dependent  problems.  We  have  shown  that  the  application  of  the  simple 
scaling  limiter,  as  is  usually  used  in  the  explicit  time  marching  case,  would  kill  the  or¬ 
der  of  accuracy  beyond  second  order,  and  we  have  designed  a  rotational  limiter  which 
can  be  shown  to  maintain  both  the  high  order  accuracy  and  positivity.  The  approach 
in  [16]  holds  a  good  potential  for  leading  to  an  efficient  and  robust  high  order  DG 
solver  for  radiative  transfer  equations.  Further  development  of  this  approach  is  now 
underway. 


2.  Broader  impact 

During  the  period  of  this  AFOSR  supported  project,  the  co-PI,  Professor  Shu, 
taught  several  undergraduate  and  graduate  level  scientific  computing  courses  covering 
related  areas,  and  directed  10  Ph.D.  students  at  Brown  University  (4  of  them  have 
already  obtained  their  Ph.D.  degrees,  2  in  2013,  1  in  2014  and  and  1  in  2015).  Among 
these  10  Ph.D.  students,  4  are  women  (women  are  underrepresented  in  mathematical 
sciences).  The  AFOSR  grant  greatly  facilitated  the  offering  of  the  courses  and  the 
training  of  these  graduate  students.  It  has  also  directly  supported  one  graduate 
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student  during  this  period. 
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1.  Background 

The  classical  finite  element  method,  based  on  low  order  piecewise  poly¬ 
nomial  approximation  using  Lagrange  shape  functions  in  conjunction  with 
mesh  refinement,  offers  geometric  flexibility  but  suffers  from  relatively  poor 
resolution.  The  spectral  method,  on  the  other  hand,  achieves  high  resolution 
by  employing  a  fixed  mesh  (often  a  single  element)  in  conjunction  with  very 
high  degree  polynomial  approximation,  but  suffers  from  a  lack  of  geomet¬ 
ric  flexibility  compared  with  the  finite  element  method.  The  spectral/fip- 
version  finite  element  method  aims  to  combine  the  advantages  of  each  ap¬ 
proach  by  using  high  degree  piecewise  polynomials  in  a  finite  element  setting. 

High  order  finite  element  methods  have  been  analysed  extensively  for  a 
wide  variety  of  applications  and  are  known  to  be  capable  of  producing  expo¬ 
nential  rates  of  convergence,  even  for  challenging  problems  with  singularities, 
sharp  boundary  layers,  and  high  frequency  oscillations.  High  order  polyno¬ 
mial  approximations  are  commonplace  in  many  areas  of  scientific  computing 
including  computer  graphics,  computer  aided-geometric  design,  and  spectral 
methods  for  PDEs.  It  is  commonplace  to  see  the  spectral  method  used  with 
approximation  orders  in  the  100s  or  1000s. 

Yet,  despite  theory  giving  the  nod  to  the  use  of  very  high  order  finite 
element  methods,  the  range  of  polynomial  degree  used  in  practical  finite  el¬ 
ement  computations  is  rarely  larger  than  eighth  order!  Possible  explanations 
for  the  use  of  comparatively  modest  polynomial  orders  include  issues  of  ef¬ 
ficiency,  implementation,  and  stability.  Moreover,  existing  implementations 
of  high  order  finite  elements  tend  to  be  memory  hungry,  often  relying  on 
the  use  of  precomputed  arrays  and  look-up  tables,  a  feature  that  does  not 
augur  well  for  the  future  given  the  nature  of  emerging  computer  hardware 
systems.  Whatever  the  underlying  reasons,  it  is  clear  that  the  rather  mod¬ 
est  polynomial  degrees  seen  in  high  order  finite  element  analysis  are  due  to 
practical  considerations  rather  than  any  theoretical  barriers. 

The  proposed  programme  of  research  is  directed  towards  tackling  these 
issues. 
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2.  Algorithms  for  Bernstein-Bezier  Finite  Elements  of  High, 
Non-Uniform  Order  in  Any  Dimension 

The  archetypal  pyramid  algorithm  is  the  de  Casteljau  algorithm;  a  stan¬ 
dard  tool  for  the  evaluation  of  Bezier  curves  and  surfaces,  and  polynomials 
expressed  in  Bernstein-Bezier  form  in  general.  The  algorithm  enjoys  wide¬ 
spread  usage,  despite  being  sub-optimal,  thanks  to  the  ease  with  which  it 
can  be  implemented,  its  underlying  stability  properties  and  short,  simple  re¬ 
cursive  nature  with  minimal  memory  overhead.  In  essence,  the  de  Casteljau 
algorithm  replaces  a  single  high  order  polynomial  by  a  recursive  sequence  of 
self-similar  affine  combinations.  Pyramid  algorithms  are  ubiquitous  in  the 
computer  aided  geometric  design  community  for  computations  using  high 
order  curves  and  surfaces.  Pyramid  algorithms  have  received  no  attention 
whatsoever  from  the  high  (or  low)  order  finite  element  community.  In  the 
article  [1]  we  develop  and  analyse  pyramid  algorithms  for  the  efficient  han¬ 
dling  of  all  of  the  basic  finite  element  building  blocks,  including  the  assembly 
of  the  element  load  vectors  and  element  stiffness  matrices. 

Hierarchic  bases  have  been  prevalent  throughout  the  high  order  finite  el¬ 
ement  literature  from  the  very  outset.  The  first  step  towards  developing 
pyramid  algorithms  for  high  order  finite  elements  is  to  dispense  with  hier¬ 
archic  bases  in  favour  of  the  non-hier archie,  Bernstein-Bezier  (BB-)  basis. 
Some  practitioners  may  baulk  at  this  prospect,  pointing  to  the  freedom 
hierarchic  bases  bring  to  allow  varying  local  approximation  order  of  the 
elements.  However,  we  shall  see  that  pyramid  algorithms  bring  the  same 
flexibility  to  the  (non-hierarchic)  BB-basis  at  no  additional  complexity ,  but, 
in  addition,  with  the  usual  advantages  associated  with  pyramid  algorithms. 
Bernstein-Bezier  bases  have  recently  been  shown  to  offer  some  advantages 
for  uniform  order  finite  element  approximation. 

We  begin  by  generalising  the  Bernstein  Decomposition  for  uniform  ap¬ 
proximation  order,  to  the  variable  order  setting.  A  non-uniform  order 
Bernstein-Bezier  basis  emerges  for  which  the  enforcement  of  conforming 
between  elements  of  differing  local  orders  is  a  natural  consequence  of  the 
structure  of  the  non-uniform  order  Bernstein  Decomposition. 

A  new,  non-uniform  order,  variant  of  the  de  Casteljau  algorithm  is  devel¬ 
oped  (expressed  as  a  pyramid  algorithm)  that  retains  the  favourable  features 
of  the  original  algorithm.  The  new  variant  is  applicable  to  the  variable  poly¬ 
nomial  order  case  but  incurs  no  additional  complexity  compared  with  the 
original  algorithm.  The  classical  degree  raising  algorithm  is  given  a  similar 
treatment  giving  a  new,  non-uniform  order,  degree  raising  pyramid  algo¬ 
rithm.  Yet  more  interestingly,  the  dual  pyramid  of  our  non-uniform  order 
degree  raising  algorithm  gives  a  pyramid  algorithm  for  the  assembly  of  the 
non-uniform  order  element  load  vector.  The  complexity  of  the  algorithm  is 
the  same  as  that  of  the  most  efficient  hierarchic  bases  currently  in  use. 

An  algorithm  for  the  construction  of  the  element  matrices  in  optimal 
complexity  for  uniform  order  approximation  on  a  simplicial  elements  was 
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developed  only  recently  [2] .  The  efficiency  of  that  algorithm  uses  the  under¬ 
lying  uniformity  of  the  polynomial  order  across  the  element  in  an  essential 
way.  There  is  no  existing  algorithm  that  constructs  the  non-uniform  order 
element  matrices  in  optimal  complexity.  In  the  final  part  of  this  work,  we 
extend  the  algorithm  of  [2]  to  the  non-uniform  order  case.  Moreover,  we 
show  that  the  resulting  algorithm  achieves  the  optimal  complexity. 

In  summary,  the  article  [1]  provides  the  methodology  that  enables  the 
efficient  use  of  a  completely  general  distribution  of  polynomial  degrees  with¬ 
out  any  restriction  in  changes  between  adjacent  cells  on  simpicial  partitions 
in  any  number  of  spatial  dimensions. 

3.  Bernstein-Bezier  FEM  on  Tet-Hex-Pyramidal  Meshes 

In  theory  at  least,  finite  element  approximation  on  a  three  dimensional 
domain  differs  little  from  approximation  over  planar  domains.  However,  the 
reality  is  that  one  is  beset  by  various  niggling  issues  not  least  of  which  is  the 
choice  of  a  partitioning  for  the  domain.  On  the  one  hand,  hexahedral  ele¬ 
ments  are  often  preferred  for  reasons  of  efficiency  stemming  from  the  tensor- 
product  nature  of  the  elements  and  the  underlying  approximation  space,  but 
the  practical  issues  of  meshing  a  complicated  three  dimensional  geometry 
using  only  hexahedra  should  not  be  underestimated.  On  the  other  hand, 
while  there  are  many  efficient  algorithms  available  for  meshing  using  tetra¬ 
hedral  elements,  the  numbers  of  degrees  of  freedom  in  the  resulting  finite 
element  space  generally  significantly  exceed  what  would  be  needed  should  a 
hexahedral  partitioning  be  available,  without  a  commensurate  improvement 
in  accuracy. 

An  obvious  solution  would  seem  to  consist  of  using  a  partition  comprised 
of  a  mix  of  both  tetrahedra  and  hexahedra  with  hexahedra  used  as  much 
as  possible,  with  the  tetrahedra  used  to  mesh  the  remaining  more  com¬ 
plicated  geometrical  features.  This  would  be  fine  were  it  not  for  the  fact 
that  tetrahedra  and  hexahedra  fail  to  tessellate.  Again,  a  possible  solution 
readily  suggests  itself  whereby  pyramidal  elements  are  used  as  a  means  of 
interfacing  between  the  tetrahedra  with  the  hexahedra. 

The  question  then  becomes  one  of  how  to  define  shape  functions  on  the 
pyramids  in  such  a  way  as  to  obtain  globally  C°  (conformal)  test  functions. 
The  problem  now  is  that  it  is  impossible  to  do  this  using  polynomial  basis 
functions.  Various  approaches  have  been  taken  to  solving  this  problem.  The 
upshot  is  that  non-polynomial  basis  functions  are  required  which,  in  turn, 
gives  rise  to  further  difficulties  associated  with  numerical  integration  and 
the  analysis. 

At  this  point,  one  might  be  left  with  a  growing  feeling  that  hybrid  parti¬ 
tions  may  be  more  trouble  than  they  are  worth.  Nevertheless,  the  current 
work  presents  a  fresh  approach  to  dealing  with  the  pyramidal  transition  ele¬ 
ments  which:  avoids  the  need  for  non-polynomial  functions  in  favour  of  using 
piecewise  polynomials,  and  as  such  is  entirely  in  keeping  with  the  spirit  of 
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the  finite  element  method;  neatly  side  steps  technical  questions  associated 
with  the  error  of  the  quadrature  and  approximation  power  of  the  spaces; 
and,  opens  up  the  possibility  for  adopting  fast  (optimal  complexity)  algo¬ 
rithms  developed  in  our  previous  work  [2]  for  the  system  matrix  assembly 
and  matrix-free  implementations. 

In  the  article  [4]  we  characterise  the  hybrid  partitions  of  interest.  While 
this  might  seem  to  be  straightforward,  in  practice,  it  is  of  considerable  inter¬ 
est  to  allow  partitions  in  which  the  faces  of  the  hexahedra  are  allowed  to  be 
bilinear  quadrilaterals.  In  turn,  this  means  that  one  is  obliged  to  consider 
pyramids  whose  quadrilateral  face  is  bilinear,  which  is  generally  disallowed 
in  existing  work  on  pyramidal  elements.  We  therefore  devote  attention  to 
carefully  defining  pyramidal  elements  of  this  type  and  characterizing  condi¬ 
tions  under  which  associated  maps  to  the  reference  element  are  invertible. 
The  bilinear  face  of  the  pyramid  means  that  the  splitting  of  the  pyramid 
into  interface  tetrahedra  is  not  straightforward.  We  discuss  the  appropriate 
way  to  split  pyramids  into  a  pair  of  interface  tetrahedra  which,  of  necessity, 
must  be  curvilinear.  We  introduce  the  approximation  space  to  be  used  in 
the  finite  element  method,  and  give  a  concrete  construction  of  the  shape 
functions  associated  with  each  of  our  elements.  We  discuss  domain  points 
and  minimal  determining  sets,  and  we  use  these  concepts  to  find  a  formula 
for  the  dimension  of  our  spline  space  and  to  construct  a  locally  supported 
basis  for  it.  Finally,  we  collect  some  useful  computational  algorithms  for 
dealing  with  Bernstein  basis  polynomials  from  our  previous  work,  which  are 
then  used  in  to  develop  algorithms  for  efficiently  computing  all  of  the  various 
quantities  needed  for  a  finite  element  analysis  and  evaluation  of  the  result¬ 
ing  approximation.  Finally,  the  approximation  power  of  our  spline  space  is 
proved. 

4.  Bernstein-Bezier  Basis  for  Arbitrary  Order  Mixed  FEM 

Bernstein  polynomials  play  an  important  role  in  CAGD  and  possess  a 
number  of  interesting  properties  that  have  led  to  their  widespread  usage 
in  a  range  of  applications.  Historically,  hierarchic  polynomial  bases  have 
been  used  virtually  exclusively  for  high  order  finite  element  approximation. 
However,  the  Bernstein  basis  is  more  natural  with  the  degrees  of  freedom 
being  more  akin  to  the  nodal  degrees  of  freedom  described  in  most  finite 
element  textbooks.  More  importantly,  it  was  recently  shown  that  using  a 
Bernstein  basis  for  the  standard  ^-conforming  finite  element  spaces  enables 
the  mass  and  stiffness  matrices  to  be  constructed  in  optimal  complexity  [2]  in 
terms  of  the  polynomial  degree.  The  algorithms  take  advantage  of  properties 
of  Bernstein  polynomials  and  the  sum- factorisation  approach  that  provides 
the  cornerstone  of  the  spectral  method. 

In  [3]  we  extend  the  ideas  of  [2]  to  the  case  of  conforming  finite  element 
approximation  of  the  space  H{ div)  on  triangulations.  We  seek  a  basis  that 
gives  a  clear  separation  between  what  are  sometimes  termed  the  curls  of 
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the  Bernstein  polynomial  basis  for  the  H1  space,  and  the  non-curls  that 
characterize  the  specific  H( div)  finite  element  space  (Raviart-Thomas  in 
our  case).  The  resulting  basis  has  two  distinct  components  reflecting  this 
separation  with  the  basis  functions  in  each  component  having  a  natural 
identification  with  a  domain  point,  or  node,  on  the  reference  element  in 
contrast  with  the  hierarchic  finite  element  bases  developed  hitherto.  More 
importantly,  it  is  shown  that  the  basis  retains  the  favourable  properties  of 
the  Bernstein  basis  that  were  used  in  [2]  to  develop  efficient  computational 
procedures  for  the  application  of  the  elements. 

5.  Solution  of  Bernstein- Vandermonde  Systems 

The  Bernstein-Bezier  (BB)  representation  of  a  polynomial  p  €  Pn([0, 1]) 
takes  the  form  p  =  Ylk=o  ck^k  which  the  coefficients  {c^}  are  referred 
to  as  control  points.  However,  while  the  polynomial  p  satisfies  p(0)  =  Co 
and  p{  1)  =  cn,  this  property  does  not  hold  at  the  remaining  control  points. 
On  the  one  hand,  this  property  does  not  hinder  the  typical  workflow  of  a 
CAGD  engineer  whereby  the  locations  of  {c^}  are  adjusted  until  a  curve 
of  the  desired  shape  is  obtained.  In  effect,  the  control  points  are  used  to 
define  the  curve  directly.  On  the  other  hand,  the  typical  usage  of  polyno¬ 
mials  in  scientific  computing  is  rather  different  in  that  one  generally  wishes 
to  fit  a  polynomial  to  a  given  function.  For  example,  in  applying  the  fi¬ 
nite  element  method  to  approximate  the  solution  of  a  partial  differential 
equation,  one  might  have  boundary  or  initial  data  and  require  to  choose  an 
appropriate  (piecewise)  polynomial  approximation  of  the  data.  The  approx¬ 
imation  is  often  chosen  to  be  an  interpolant,  leading  to  what  we  shall  term 
the  Bernstein-Bezier  interpolation  problem,  which  consists  of  computing  the 
control  points  {cfc}j!=0  such  that  the  associated  Bernstein-Bezier  polynomial 
interpolates  the  data: 


n 

P(xj)  =  °kBk(xj)  =  fji  j  =  0,  •  •  • ,  n.  (1) 

k= 0 

Conditions  (1)  may  be  expressed  as  a  system  of  linear  equations  involv¬ 
ing  the  Bernstein-  Vandermonde  matrix.  If  the  monomial  basis  were  to  be 
used,  then  the  standard  Vandermonde  matrix  would  emerge.  The  highly  ill- 
conditioned  nature  of  the  Vandermonde  matrix  is  well-documented.  Notwith¬ 
standing,  the  inversion  of  the  Vandermonde  matrix  to  compute  the  La- 
grangian  interpolant  is  in  some  ways  preferable  to  more  direct  methods. 

The  Bernstein- Vandermonde  matrix  is  likewise  found  to  be  highly  ill- 
conditioned,  though  to  a  lesser  extent  than  the  Vandermonde  matrix,  sug¬ 
gesting  that  its  inversion  may  not  be  the  wisest  approach.  Nevertheless, 
structure  of  the  matrix  arising  from  the  total  positivity  of  the  Bernstein  ba¬ 
sis  means  that  using  Neville  elimination  to  solve  the  system  obviates  some 
of  the  issues  due  to  ill-conditioning.  Marco  and  Martinez  exploit  this  fact 
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and  derive  an  algorithm  for  the  inversion  of  the  matrix  that  has  0(n2) 
complexity — the  same  as  multiplying  by  the  inverse  of  the  matrix. 

Remarkable  though  the  Marco-Martfnez  algorithm  may  be,  it  does  have 
its  drawbacks: 

•  the  derivation  of  the  algorithm  is  highly  technical  involving  non¬ 
trivial  identities  for  the  minors  of  Bernstein- Vandermonde  matrices; 

•  the  interplay  between  the  ideas  of  total  positivity,  bidiagonal  factori¬ 
sation  and  Neville  elimination  are  not  part  of  the  standard  armoury 
of  many  non-specialists; 

•  the  algorithm  seems  to  be  firmly  rooted  to  the  solution  of  the  Bern¬ 
stein  interpolation  problem  in  the  univariate  case  when  simplices 
are  considered  (indeed,  total  positivity  is  essentially  a  univariate 
concept). 

In  [5]  we  address  these  issues.  Specifically,  we  present  an  alternative 
algorithm  for  the  inversion  of  the  univariate  Bernstein- Vandermonde  system 
that  has: 

•  the  same  complexity  as  the  Marco-Martfnez  algorithm  and  whose 
stability  does  not  seem  to  be  in  any  way  inferior; 

•  a  simple  derivation  that  could  be  taught  to  undergraduates  familiar 
with  only  the  basic  theory  of  Lagrange  interpolation  (at  least  in  the 
univariate  case); 

•  a  natural  generalisation  to  the  multivariate  case  on  simplices  (essen¬ 
tial  for  applications  such  as  finite  element  analysis  in  two  and  three 
dimensions). 

More  specifically,  we  derive  the  new  algorithm  in  the  univariate  case  using 
only  elementary  facts  about  univariate  interpolation  and  basic  properties 
of  Bernstein  polynomials,  and  present  a  straightforward  extension  of  the 
univariate  algorithm  to  the  multivariate  tensor  product  case.  We  tackle 
the  much  harder  task  of  solving  the  Bernstein  interpolation  problem  on 
simplices  in  higher  dimensions  in  which  even  the  existence  of  the  Lagrangian 
interpolant  is  less  obvious.  Using  the  standard  hypothesis  under  which  the 
Lagrange  interpolant  exists,  we  develop  an  algorithm  for  the  solution  of  the 
multivariate  Bernstein  interpolation  problem  on  simplices  in  two  dimensions, 
and  indicate  how  it  may  be  extended  to  arbitrary  dimensions.  The  ideas  used 
are,  modulo  technical  issues,  essentially  the  same  to  those  used  in  Section  2 
to  handle  the  univariate  case.  We  give  numerical  examples  illustrating  the 
behaviour  and  stability  of  the  resulting  algorithms. 

The  algorithms  developed  in  the  present  work  form  a  key  step  towards  a 
more  widespread  use  of  Bernstein-Bezier  techniques  in  scientific  computing 
in  general,  and  in  finite  element  analysis  in  particular  by  enabling  the  use  of 
non-homogeneous  initial  and  boundary  data.  More  generally,  the  problem 
of  how  to  solve  Bernstein- Vandermonde  systems  on  simplices  in  arbitrary 
dimension  is  addressed. 
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1  Stability  and  Numerics  for  Weakly  Interacting 
Particle  Systems 

1.1  Problem  formulation  and  issues 

Many  physical  and  engineering  systems  can  be  modeled  as  a  large  collection 
of  stochastically  evolving  agents  or  particles,  whose  dynamics  are  weakly 
coupled  by  an  interaction  that  depends  only  on  the  distribution  of  the  par¬ 
ticles.  Simulations  of  these  systems  has  shown  that  the  dynamics  of  these 
systems  exhibit  interesting  and  complex  behavior.  The  number  of  particles 
is  typically  large,  and  there  are  no  general  purpose  tools  for  analyzing  such 
systems  that  are  both  accurate  and  computationaly  feasible.  Of  particu¬ 
lar  interest  are  the  long-term  stability  and  nretastability  properties  of  these 
systems. 

This  part  of  the  research  project  was  concerned  with  the  development  of 
general  mathematical  methods,  both  analytical  and  numerical,  to  study  the 
long-time  behavior,  nretastability  properties  and  control  of  such  stochastic 
systems.  In  many  cases,  the  N  particles  are  exchangeable  and  the  relevant 
quantities  are  captured  by  the  empirical  distribution  of  the  A-particle  sys¬ 
tem.  We  studied  two  classes  of  problems,  one  when  the  state  of  each  particle 
is  discrete,  taking  a  finite  number  of  values,  and  when  the  state  of  each  par¬ 
ticle  is  a  continuum,  for  example,  the  whole  real  line.  As  described  below, 
the  methods  used  in  each  of  the  two  cases  are  somewhat  different. 

1.2  Results  and  Key  Findings:  The  Finite-Dimensional  Case 

1.2.1  Accelerated  Monte  Carlo  and  large  and  moderate  devia¬ 
tions 

One  focus  was  on  the  use  of  large  deviation  methods  for  accelerated  Monte 
Carlo  schemes.  Specifically,  we  used  large  deviation  theory  for  purposes  of 
design  and  analysis,  and  when  existing  theory  is  not  adequate,  it  was  devel¬ 
oped  in  the  form  suited  to  such  applications.  The  large  deviation  theory  for 
various  process  models  was  developed  in  [BCD13,  DJ15]  and  [BDG16].  It  is 
known  that  when  using  any  of  the  standard  methods  that  reduce  variance 
and  thereby  improve  efficiency  one  needs  some  information  on  the  condi¬ 
tional  distributions  of  the  associated  process  model.  Also,  in  the  setting  of 
rare  events  the  use  of  accelerated  Monte  Carlo  is  essential,  since  the  relative 
variance  (i.e. ,  the  variance  of  the  sample  to  the  size  of  the  probability  being 
estimated)  scales  inversely  to  the  size  of  the  probability. 
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The  papers  [DJ15]  and  [BDG16]  prove  what  are  called  moderate  de¬ 
viation  principles  (or  moderate  deviation  approximations)  for  fairly  com¬ 
plicated  process  models.  A  moderate  deviation  principle  gives  information 
that  is  not  as  far  out  in  the  tails  of  the  distribution  as  in  the  case  of  large 
deviation,  and  hence  can  be  used  to  study  events  that  are  rare  but  not 
extremely  small.  When  using  a  moderate  deviation  approximation  one  is 
typically  interested  in  accelerated  methods  because  the  expense  associated 
with  generating  even  a  single  sample  can  be  high. 

As  an  example  of  how  the  results  of  [BDG16]  could  be  used  consider  the 
problem  of  assessing  the  probability  of  exceeding  an  allowed  pollution  level 
in  the  context  of  groundwater  contamination.  Pollutants  are  modeled  as 
entering  a  waterway  according  to  a  spatially  distributed  Poisson  input,  after 
which  they  diffuse  throughout  the  waterway.  The  problem  of  interest  might 
be  a  rare  event,  such  as  exceeding  a  regulatory  threshold  for  chlorinated 
hydrocarbons,  or  a  less  rare  event,  such  as  higher-than-normal  levels  of  plant 
nutrients,  leading  to  algae  bloom  and  fish  kill.  Each  sample  thus  requires 
the  solution  of  a  2  or  3  dimensional  partial  differential  equation. 

The  process  models  covered  in  [DJ15]  are  more  directly  relevant  to  the 
project’s  overall  goals,  and  include  weakly  interacting  many  particle  models 
with  a  finite  dimensional  empirical  measure  (e.g.,  weakly  interacting  finite 
state  Markov  chains).  Here  again  one  is  especially  concerned  with  variance 
reduction  due  to  the  fact  that  simulating  even  a  single  sample  is  costly  when 
the  number  of  particles  is  large. 

For  both  the  large  deviation  and  moderate  deviation  approaches  the  con¬ 
struction  of  accelerated  Monte  Carlo  schemes  are  based  on  subsolutions  to  a 
nonlinear  partial  differential  equation  (PDE)  associated  with  the  rate  func¬ 
tion  of  the  process.  It  is  the  rate  function  which  gives  approximations  to 
conditional  distributions  of  the  process.  In  general  one  expects  the  moderate 
deviations  rate  function  to  involve  a  linear  approximation  of  the  dynamics 
and  a  quadratic  approximation  of  the  costs  found  in  the  large  deviations  rate 
function,  centered  around  the  law  of  large  numbers  limit.  This  corresponds 
to  a  “local  Gaussian”  approximation,  and  the  PDEs  involved  are  those  as¬ 
sociated  with  the  famous  “linear-quadratic-regulator,”  whose  solution  can 
be  constructed  in  terms  of  the  solution  to  an  appropriate  Riccati  equation. 

In  [DJ16]  the  moderate  deviation  results  of  [DJ15]  are  applied  to  the 
design  of  efficient  importance  sampling  schemes.  We  present  various  imple¬ 
mentations,  and  include  numerical  results  to  contrast  their  performances, 
including  weakly  interacting  particle  models.  Also  included  is  a  discussion 
under  what  circumstances  a  particular  scaling  (large  or  moderate)  is  most 
appropriate.  Owing  to  the  fact  that  the  linear-quadratic-regulator  has  (up 
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to  the  solution  of  a  Riccati  equation)  has  an  explicit  solution,  the  design  of 
optimal  and  nearly  optimal  schemes  can  be  done  in  great  generality. 


1.2.2  Robustness  and  model  error  in  the  context  of  rare  events 

There  is  little  known  regarding  the  robustness  of  large  deviation  estimates 
with  regard  to  modeling  errors.  This  is  the  problem  treated  in  [ACD15], 
which  builds  on  ideas  first  developed  in  [CD13] .  In  [ACD15],  we  extend 
the  duality  between  exponential  integrals  and  relative  entropy  to  a  varia¬ 
tional  formula  for  exponential  integrals  involving  the  Renyi  divergence.  This 
formula  characterizes  the  dependence  of  risk-sensitive  functionals  to  pertur¬ 
bations  in  the  underlying  distribution.  It  also  shows  that  perturbations  of 
related  quantities  determined  by  tail  behavior,  such  as  probabilities  of  rare 
events,  can  be  bounded  in  terms  of  the  Renyi  divergence.  The  character¬ 
ization  gives  rise  to  tight  upper  and  lower  bounds  that  are  meaningful  for 
all  values  of  a  large  deviation  scaling  parameter,  allowing  one  to  quantify 
in  explicit  terms  the  robustness  of  risk-sensitive  costs.  As  applications  we 
consider  problems  of  uncertainty  quantification  when  aspects  of  the  model 
are  not  fully  known,  as  well  their  use  in  bounding  tail  properties  of  an  in¬ 
tractable  model  in  terms  of  a  tractable  one. 

1.2.3  Cooperative  control  of  many  particle  systems 

In  [DRV16]  we  study  many  particle  exit  time  stochastic  control  problems 
with  risk-sensitive  cost.  Each  particle  takes  values  in  a  finite  set  A,  and 
by  controlling  the  rates  of  each  individual  the  aim  is  to  keep  the  system 
away  from  a  “ruin”  set  K,  for  as  long  as  possible  and  with  the  least  cost.  We 
prove,  under  suitable  assumptions,  that  for  every  finite  number  n  of  particles 
the  control  problem  is  equivalent  to  one  with  an  ordinary  (additive)  cost. 
Moreover,  when  K,  C  Xn  can  be  identified  with  a  subset  of  the  simplex  of 
probability  measures  V(X)  (in  the  sense  that  for  every  permutation  a  E  Sn 
we  have  a  1C  =  1C),  then  we  can  replace  the  original  problem  by  one  on 
Vn(X)  =  V(X)  (~l  jTC1 ,  getting  in  this  way  a  control  problem  whose  state  is 
the  empirical  measure  on  the  states  of  the  individual  particles. 

We  start  with  a  collection  of  independent  particles/agents  with  each 
evolving  according  to  the  transition  rates  7.  This  is  the  “preferred”  or 
“nominal”  dynamic,  and  is  what  would  occur  if  no  “outside  influence”  or 
other  form  of  control  acts  on  the  particles.  If  a  controller  should  wish 
to  change  this  behavior,  then  it  must  expect  to  pay  a  cost  to  do  so.  We 
model  the  situation  where  limited  information  on  the  system  state,  and  in 
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particular  information  relating  only  to  the  empirical  measure  of  the  states  of 
all  particles,  is  used  to  produce  a  desired  behavior  of  the  group  of  particles, 
which  again  will  be  characterized  in  terms  of  their  empirical  measure  (how 
the  collective  “loads”  the  system). 

The  paper  presents  three  results.  The  first  is  that,  under  certain  con¬ 
ditions,  for  each  n  the  risk-sensitive  control  problem  is  equivalent  to  an 
ordinary  control  problem  rather  than  a  stochastic  game ,  as  one  might  ex¬ 
pect.  We  also  show  under  appropriate  conditions  that  both  the  risk  sensitive 
and  the  ordinary  control  problems  are  equivalent  to  mean  field  control  prob¬ 
lems.  The  last  contribution  is  convergence  of  a  suitably  normalized  value 
functions  for  the  n  particle  system  to  the  value  function  of  a  deterministic 
control  problem. 

1.3  Results  and  Key  Findings:  The  Infinite-dimensional  Case 
1.3.1  Network  models 

Large-scale  stochastic  networks  arise  in  a  variety  of  of  real  world  appli¬ 
cations,  ranging  from  telecommunications,  service  systems  and  computer 
networks  to  health  care  services  and  biological  systems.  Such  networks 
are  typically  too  complex  and  not  amenable  to  exact  analysis.  However 
it  is  often  possible  to  provide  insight  into  network  performance,  in  both 
the  transient  and  equilibrium  regimes,  through  tractable  approximations. 
Finite-dimensional  Markov  processes  have  been  used  extensively  in  the  liter¬ 
ature  to  describe  the  scaling  limits  of  stochastic  networks  under  simplifiying 
assumptions  such  as  exponential  service  distributions.  However,  statistical 
analysis  has  shown  that  the  service  distribution  is  typically  non-exponential. 

In  the  works  [A15]  and  [A16a],  the  authors  study  the  classical  model 
of  an  N- server  queue  with  cumulative  renewal  arrivals  and  jobs  with  i.i.d. 
service  distribution  served  in  a  First-Come-First-Serve  (FCFS)  manner.  It 
was  shown  in  earlier  work  of  Kaspi  and  Ramanan  (2011)  that  such  a  system 
can  be  viewed  as  a  weakly  interacting  particle  system  in  which  each  particle 
corresponds  to  a  server,  and  the  state  of  the  system  is  the  amount  of  time 
that  the  server  has  spent  serving  the  job  it  currently  has  in  system  (or,  is  in 
an  idle  state,  if  it  is  not  serving  any  job).  Since  the  servers  are  exchangeable, 
the  dynamics  is  captured  by  the  empirical  measure  of  the  state  of  the  servers. 
Kaspi  and  Ramanan  (2011)  established  a  certain  functional  laws  of  numbers 
limit  for  this  measure-valued  system,  and  subsequently,  in  2013,  established 
a  functional  central  limit  theorem.  Both  of  these  limit  theorems  captured 
the  mean  behavior  and  fluctuations  around  the  mean  behavior  on  finite  time 
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intervals.  A  major  open  problem  posed  by  Halfin  and  Whitt  in  1981,  was 
to  obtain  an  approximation  for  the  diffusion-scaled  steady  state  distribution 
of  this  network.  A  key  challenge  was  that,  due  to  the  infinite-dimensional 
of  the  state-space,  classical  Lyapunov  function  and  positive  Harris  recur¬ 
rence  techniques  are  not  applicable.  In  [A15]  and  [A16a]  we  introduce  new 
methods  that  rely  on  the  framework  of  asymptotic  coupling  recently  intro¬ 
duced  by  Hairer,  Mattingly  and  Scheutzow  and  others  to  show  (for  a  large 
class  of  service  distributions)  convergence  of  the  diffusion  scaled  N- server 
queue  lengths  to  a  limit,  that  can  be  characterized  as  the  unique  stationary 
distribution  of  an  infinite-dimensional  Markov  process. 

The  methods  introduced  in  [A15]  and  [A16a],  in  particular  certain  asymp¬ 
totic  coupling  constructions,  are  more  broadly  applicable  for  study  of  a  much 
broader  class  of  networks,  and  provide  a  major  first  step  in  developing  an¬ 
alytical  characterizations  and  numerical  methods  for  estimating  the  limit 
steady-state  distribution. 

The  work  in  [A16b]  studies  a  more  complicated  many-server  model  with 
randomized  load  balancing,  in  which  each  server  has  its  own  dedicated 
queue.  This  system  can  once  again  be  represented  by  a  weakly  interacting 
particle  system,  whose  dyanamics  is  captured  by  interacting  measure- valued 
processes.  This  randomized  load  balancing  network  had  earlier  only  been 
analyzed  in  the  case  of  exponential  service  distributions  independently  by 
Mitzenmacher  and  by  Dobrushin,  Karpelevich  and  Vvydenskaya.  Our  work 
provides  a  framework  for  analyzing  such  models  in  the  presence  of  general 
service  distributions.  In  particular,  in  [A16b]  we  show  that,  as  the  num¬ 
ber  of  servers  N  goes  to  infinity,  the  measure-valued  process  that  captures 
the  dynamics  of  the  load  balancing  network,  converges  to  a  hydrodynamic 
limit.  This  allows  one  to  study  the  impact  of  the  service  distribution  on  the 
performance.  Our  general  framework  is  also  useful  for  the  study  of  a  much 
broader  class  of  many-server  load  balancing  networks. 

1.3.2  Particles  with  Singular  Interaction 

There  is  a  substantial  interest  in  interacting  particle  systems  in  which  each 
particle  takes  values  in  the  reals,  and  the  n-particle  configurations  follow  a 
Gibbs  distribution,  which  is  parameterized  by  a  quantity  (often  referred  to 
as  the  inverse  temperature  (3n)  and  an  energy  functional  that  is  the  sum  of 
a  (possibly  singular)  interaction  W  and  a  confining  potential  V .  Such  pro¬ 
cesses  model  a  variety  of  systems  of  interest  including  interacting  Brownian 
particles,  eigenvalues  of  random  matrices,  and  arise  in  the  study  of  sampling 
and  simulated  annealing.  In  particular,  it  is  of  interest  to  study  probabilties 
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of  rare  events  in  such  systems,  when  the  number  of  particles  is  large.  In 
[DLR15]  we  establish  large  deviation  principles  (LDPs)  for  empirical  mea¬ 
sures  associated  with  a  sequence  of  Gibbs  distributions  of  n-particle  config¬ 
urations  associated  with  the  parameter  / 3n ,  both  with  speeds  (3n/n  — >  oo, 

While  the  case  when  the  interaction  W  is  bounded  is  well-studied,  in 
the  applications  mentioned  above,  W  may  blow  up  along  the  diagonal.  In 
[DLR15],  in  which  case  the  rate  function  is  expressed  in  terms  of  a  func¬ 
tional  involving  the  potentials,  and  with  the  speed  f3n  =  n,  when  the  rate 
function  contains  an  additional  entropic  term.  Using  the  weak  convergence 
methods  first  introduced  by  Dupuis  and  Ellis  we  establish  large  deviation 
principles  not  only  with  respect  to  the  weak  topology,  but  also  with  respect 
to  stronger,  Wasserstein-type  topologies,  which  allows  one  to  also  under¬ 
stand  the  behavior  of  unbounded  continuous  functionals  of  the  empirical 
measure  of  the  particles,  such  as  moments,  which  are  of  interest  in  appli¬ 
cations.  Our  work  provides  a  common  framework  for  the  analysis  of  LDPs 
with  all  speeds,  and  includes  cases  not  covered  due  to  technical  reasons  in 
previous  work  of  Ben  Arous  et.  al  (1998)  and  Chafai  (2014). 

1.3.3  Sensitivities  of  Constrained  Processes 

Reflected  stochastic  processes  that  are  constrained  to  lie  in  the  closure  of 
a  convex  polyhedral  domain  arise  in  many  contexts,  including  as  diffusion 
approximations  of  stochastic  networks,  in  the  study  of  interacting  diffusions 
and  limits  of  interacting  particle  systems,  in  chemical  and  biochemical  re¬ 
action  networks  and  in  mathematical  finance.  The  analysis  of  processes 
with  state  constraints  is  challenging  due  to  the  fact  that  the  nature  of  the 
dynamics  often  changes  abruptly  at  the  boundary  of  the  domain,  and  is 
further  complicated  when  the  boundary  is  not  smooth.  Given  that  in  many 
models,  there  is  uncertainty  in  the  parameters  that  define  the  system,  a  key 
question  of  interest  is  the  dependence  of  the  reflected  process  (in  particular, 
expectations  of  functionals  of  the  process)  on  various  parameters  that  define 
it. 

The  sensitivity  of  the  stochastic  process  to  perturbations  in  the  initial 
condition  and  other  parameters  that  define  the  process,  is  a  classical  topic 
in  stochastic  analysis.  For  example,  there  is  a  substantial  body  of  work  that 
studies  these  questions  for  (unconstrained)  diffusions  in  Euclidean  space, 
with  contributions  from  Elworthy  (1978),  Bismut  (1981),  Ikeda  and  Watan- 
abe  (1981),  Kunita  (1981),  Metivier  (1982)  and  others.  In  contrast,  there 
are  relatively  few  results  for  reflected  processes  or  even  reflected  Brownian 
motions,  especially  in  the  context  of  oblique  reflection  and  nonsmooth  do- 
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mains  which  is  relevant  in  applications.  A  convenient  tool  for  the  analysis  of 
reflected  Brownian  motions  is  the  so-called  extended  Skorokhod  map  (ESM), 
which  maps  an  unconstrained  trajectory  to  a  version  that  is  constrained  to 
lie  in  a  certain  specified  domain  by  applying  a  certain  constraining  force 
in  the  (possibly)  oblique  reflection  vector  field  specified  on  the  boundary. 
In  [LR15]  we  first  show  that  the  analysis  of  pathwise  differentiability  of  re¬ 
flected  Brownian  motions  can  be  reduced  to  the  study  of  so-called  directional 
derivatives  of  the  associated  ESM.  For  a  special  class  of  ESMs  associated 
with  the  orthant  and  inwardly  pointing  directions  of  reflection,  which  arise 
in  many  applications,  directional  derivatives  of  the  ESM  were  studied  in 
previous  work  of  Mandelbaum  and  Ramanan  (2010).  This  class  of  ESMs 
are  particular  tractable  because  they  exhibit  a  certain  monotonicity  prop¬ 
erty  and  admit  a  semi-explicit  characterization.  However,  for  more  general 
ESMs  that  capture  a  broader  range  of  applications,  this  monotonicity  prop¬ 
erty  does  not  hold  and  the  approach  of  Mandelbaum  and  Ramanan  can  no 
longer  be  applied.  In  [LR15]  we  develop  a  completely  new  axiomatic  frame¬ 
work  for  the  characterization  of  directional  derivatives  of  a  much  broader 
class  of  ESMs  that  are  Lipschitz  continuous.  This  work  is  a  key  step  in 
future  work  that  will  study  sensitivity  analysis  of  functionals  of  RBMs,  over 
both  the  finite  and  infinite  time  horizon. 
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ACD  Atar,  R.,  Chowdhary,  K.  and  Dupuis,  P.  Robust  bounds  on  risk- 
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25,  2015,  2909-2958. 
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maps  in  convex  polyhedral  domains,  Preprint,  http://arxiv.org/ 
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A16  R.  Aghajani.  Infinite-dimensional  scaling  limits  of  stochastic  networks, 
PhD  Thesis,  Brown  University,  May  2016. 
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a  many-server  queue,  Preprint,  http://arxiv.org/abs/1512.02929 
2015. 

AR16a  R.  Aghajani  and  Ramanan,  K.  The  limit  of  stationary  distributions 
of  many-server  queues  in  the  Halfin- Whitt  regime  Preprint,  https : 
//arxiv.org/abs/1610. 01118  2016. 
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AR16b  R.  Aghajani  and  Ramanan,  K.  Hydrodynamic  limits  of  randomized 
load-balancing  networks  Preprint,  2016. 


BDG16  Budhiraja,  A.,  Ganguly,  A.  and  Dupuis,  P.  Moderate  deviation  princi¬ 
ples  for  stochastic  differential  equations  with  jumps,  to  appear  in  the 
Annals  of  Probability. 

DJ16  Dupuis,  P.  and  Johnson,  D.  Moderate  deviations  based  importance 
sampling  for  stochastic  recursive  algorithms,  to  appear  in  the  J.  of 
Applied  Probability. 

DLR16  Dupuis,  P.,  Laschos,  V.  and  Ramanan,  K.  Risk-sensitive  cooperative 
control  for  a  many  particle  system,  In  preparation,  2016. 


10 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


AI-OSR  Deliverables  Submission  Survey 


Response  ID:7089  Data 

1. 

Report  Type 
Final  Report 
Primary  Contact  Email 

Contact  email  if  there  is  a  problem  with  the  report. 

dupuis@dam.brown.edu 

Primary  Contact  Phone  Number 

Contact  phone  number  if  there  is  a  problem  with  the  report 

401  863  3238 

Organization  /  Institution  name 

Brown  University 

Grant/Contract  Title 

The  full  title  of  the  funded  effort. 

Methods  for  high-order  multi-scale  and  stochastic  problems:  analysis,  algorithms  and  applications 

Grant/Contract  Number 

AFOSR  assigned  control  number.  It  must  begin  with  "FA9550"  or  "F49620"  or  "FA2386". 

FA9550-1 2-1 -0399 

Principal  Investigator  Name 

The  full  name  of  the  principal  investigator  on  the  grant  or  contract. 

Paul  Dupuis 
Program  Officer 

The  AFOSR  Program  Officer  currently  assigned  to  the  award 
Jean-Luc  Cambier 
Reporting  Period  Start  Date 

07/1 5/2012 

Reporting  Period  End  Date 

07/1 4/2016 

Abstract 

The  project  focused  on  the  development,  analysis,  implementation,  and  application  of  efficient  and  high- 
order  accurate  methods  for 

multi-scale  and  stochastic  problems.  Research  focused  on  three  topics: 

(1 )  High  order  weighted  essentially  non-oscillatory  finite  difference  and  finite  volume  schemes, 
discontinuous  Galerkin  finite  element 

method,  and  related  methods,  for  solving  computational  fluid  dynamics  (CFD)  problems  and  other 
applications  containing  strong  shock  waves  and  complicated  smooth  region  structures. 

(2)  Development  of  higher  order  piecewise  polynomial  approximation  for  finite  element  methods. 

(3)  The  development  of  methods  of  simulation  and  analysis  for  the  study  of  large  scale  stochastic  systems 
of  weakly  interacting  particles. 

Distribution  Statement 

This  is  block  12  on  the  SF298  form. 

Distribution  A  -  A I tSi p roved  for  public  release. 


Explanation  for  Distribution  Statement 

If  this  is  not  approved  for  public  release,  please  provide  a  short  explanation.  E.g.,  contains  proprietary  information. 
SF298  Form 

Please  attach  your  SF298  form.  A  blank  SF298  can  be  found  here.  Please  do  not  password  protect  or  secure  the  PDF 
The  maximum  file  size  for  an  SF298  is  50MB. 

PD_AFOSR_sf0298.pdf 

Upload  the  Report  Document.  File  must  be  a  PDF.  Please  do  not  password  protect  or  secure  the  PDF  .  The 
maximum  file  size  for  the  Report  Document  is  50MB. 

combined.pdf 

Upload  a  Report  Document,  if  any.  The  maximum  file  size  for  the  Report  Document  is  50MB. 

Archival  Publications  (published)  during  reporting  period: 

Please  see  the  report 

New  discoveries,  inventions,  or  patent  disclosures: 

Do  you  have  any  discoveries,  inventions,  or  patent  disclosures  to  report  for  this  period? 

No 

Please  describe  and  include  any  notable  dates 

Do  you  plan  to  pursue  a  claim  for  personal  or  organizational  intellectual  property? 

Changes  in  research  objectives  (if  any): 

Change  in  AFOSR  Program  Officer,  if  any: 

The  program  officer  was  changed  from  Fariba  Fahroo  to  Jean-Luc  Cambier 
Extensions  granted  or  milestones  slipped,  if  any: 

AFOSR  LRIR  Number 
LRIR  Title 
Reporting  Period 
Laboratory  Task  Manager 
Program  Officer 
Research  Objectives 


Technical  Summary 

Funding  Summary  by  Cost  Category  (by  FY,  $K) 


Starting  FY 

FY+1 

FY+2 

Salary 

Equipment/Facilities 

Supplies 

Total 

Report  Document 

Report  Document  -  Text  Analysis 

Report  Document  -  Text  Analysis 

Appendix  Documents 

2.  Thank  You 

E-mail  user 

Oct  12,  2016  08:41 :16  Success:  Email  Sent  to:  dupuis@dam.brown.edu 
DISTRIBUTION  A:  Distribution  approved  for  public  release. 


